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ABSTRACT 


A finite element shallow-water model is tested with two types of 
surface topography. The model uses rectangular subdivisions in a vorticity- 
divergence formulation, and a semi-implicit time discretization. In the first 
experiment an east-west ridge or valley is placed in a channel with east- 
west periodic conditions. Linear quasi-geostrophic solutions are derived 
with the rigid lid assumption. The Rossby waves are successfully simulated 
in the model with linear solutions as the initial conditions. The model phase 
speeds are very close to the analytic values when the latter are properly 
corrected. In the second experiment a ridge is placed across the channel and 
the Coriolis parameter is set to zero. The initial conditions consist of a 
uniform flow through the channel and constant free-surface height. The 
numerical simulations agree with hydraulic jump theory. In the jump cases 
the model predicts increasing wind speeds and decreasing free surface 
heights. Higher spatial resolution would be required to properly simulate 


the details of the hydraulic jump formation. 
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I. INTRODUCTION 


A. GENERAL 

Wilhelm Bjerkens (1904) was the first to point out that the future 
meteorological conditions can in principle be obtained by an integration of 
differential equations which govern the behavior of the atmosphere. Such 


an integration performed using numerical methods is called numerical 


Richardson was the first to attempt a numerical weather prediction. 
After very long and time-consuming computations, he obtained a totally 
unacceptable result (Richardson, 1922). That wrong result, and 
Richardson's estimate that 64,000 men are required to advance the 
calculations as fast as the weather itself 1s advancing, left some doubt 
about the practical use of the method. However, a number of developments 
that followed improved the situation. Mainly due to the work of Rossby in 
the late 1930's, it became clear that even a rather simple equation, one that 
describes the conservation of absolute vorticity following the motion of air 
particles, suffices for an approximate description of large scale motions of 
the atmosphere. 

Finally, in 1945, the first electronic computer, ENIAC, was 
constructed. The absolute vorticity conservation equation, and ENIAC, 
were used by Charney, Fjortoft and Von Neumann in the late 1940's for the 


first successful numerical forecast (Charney et al. 1950). Much faster 


computers, and improved understanding of computational problems, now 
also enable long-term integrations of the basic primitive equations. With 
the introduction of each new generation of computers, the gap between 
numerical forecasts and atmospheric observations has decreased, but the 
rate at which this gap is decreasing appears to be leveling off. This 
indicates that technological improvements in computing power may not be 
the primary limitation to better numerical forecasts. 

During the past 15 years, there has been a significant effort within the 
numerical weather prediction area in developing limited-area, fine-mesh 
primitive equations models and applying them to operational, short range 
weather forecasts. An important practical motivation for the development of 
regional models has been the limited success of operational global models 
in the prediction of precipitation and severe weather. A parallel motivation 
for the development of regional models is their potential scientific value to 
researchers studying the structure and dynamics of mesoscale phenomena. 


(Keyser and Uccellini, 1987) 


B. BACKGROUND 

There are two fundamental methods of simulating the atmospheric 
flow; physical-models and mathematical-models techniques. With the 
physical-models technique, we construct scale model replicas of observed 
ground surface characteristics and insert them into a chamber such as a 


wind tunnel. The flow of air in this chamber is adjusted so as to best 


represent the large scale, observed atmospheric conditions. Mathematical 
modeling, on the other hand, makes use of such basic analysis techniques 
as algebra and calculus to solve directly the equations governing 
atmospheric flows. 

Numerical solution of the equations of motion is performed using the 
grid point method for most applications. Following this method a set of 
points is introduced in the region of interest and dependent variables are 
initially defined and subsequently computed at these points. This set of 
points is called a grid. It is very important to note that most of the time the 
grid points are at fixed locations in the horizontal. This means that, 
according to the Eulerian system of equations, space and time coordinates 
are chosen as independent variables. 

With the grid point method, the most common way of solving the 
governing equations is to find approximate expressions for derivatives 
appearing in the equations. The required approximate expressions are 
defined using values of the dependent variables only defined at the grid 
points, and at discrete time intervals. Thus, they are formed using 
differences of dependent variables over finite space and time intervals ; that 
is the reason why this approach is called the finite difference method. The 
approximations for derivatives are then used to construct a system of 
algebraic equations that approximates the governing partial differential 
equations. This set of algebraic equations is to be solved, usually using an 


electronic computer, by a proper step-wise procedure in time. 


A major limiting factor of finite difference approximations is the 
truncation error. That is, the smaller the grid interval, the smaller the 
truncation error. For a finite difference model to improve its accuracy, it 
would require increasing the grid matrix density. This would require 
increased computer storage and computational time. The problem here is 
not a simple one. It goes beyond numerical techniques and computer 
technology. For instance if we further reduce the grid spacing on the 7LPE 
(National Weather Service 7 Layer Primitive Equation Model) we do not get 
any significant improvement in the accuracy of the solution (Woodward, 
1981). The required additional computer capability can not be utilized using 
finite difference methods. Therefore, new numerical integration techniques 
must be investigated. 

Two alternative techniques, the spectral method and the finite element 
method, have been the subject of intensive research. Both the spectral and 
finite element methods require more computational time per forecast time 
step than does the finite difference method. That is why both the spectral 
and finite element methods must utilize efficient numerical techniques to be 
considered a viable option for numerical weather prediction. For long range 
weather predictions, the spectral method appears to be a natural method, 
when applied over the globe or hemisphere. This is due to the existence of 
efficient transforms for the nonlinear terms in spherical geometry. 


However, the spherical harmonics are globally rather than locally defined. 


That is the reason why the finite element method appears to be more 
suitable for problems of more detailed limited area forecasting. 

The Galerkin Finite Element Method (GFEM) has the potential to 
increase efficiently the spatial resolution for the purpose of simulating 
accurately the small-scale processes. More specifically, the GFEM model 
used by Hinsman, 1983, has demonstrated desirable characteristics. It 
propagates atmospheric waves better than an equivalent finite difference 
model. It also allows variable-resolution grids and responds better than an 
equivalent finite difference model near the smallest gridlength. Moving 
grids can be achieved with no apparent noise generation. Finally, it can 
utilize direct solvers and is a natural choice for vectorization on large 
computers. Furthermore, the superior small-scale response of the GFEM 
indicates potential increase in skill for regional forecasting, and it truly is a 


viable option for simulation of atmospheric flow. 


C. OBJECTIVES 

The main objective of this study is to test the Galerkin Finite Element 
model which was developed by Hinsman (1983), with topography. This 
shallow-water model uses bilinear basis functions. An earlier study with 
the same model was carried out by Neta et al. (1986), in which the bottom 
topography changes linearly to the north. 

In this thesis two types of bottom topography will be considered in a 


channel domain encompassed by solid north-south walls with east-west 


boundary conditions. In the first case the bottom is composed of two 
regions. These regions have a constant but opposite northward slope so 
that they can form either an east-west oriented ridge or an east-west 
oriented valley in the surface topography. Lines of constant bottom height 
run parallel to the x-axis, and pure geostrophic motion 1s possible only if 
the v-component of the velocity is identically zero. If we consider no mean 
flow and no beta effect, topographic Rossby waves can exist moving either 
east or west depending on the slope of the bottom. 

In the second case, a mean zonal flow passes over a ridge which 
extends north-south across the channel. It is clear now that the lines of 
constant bottom height run parallel to the y-axis. The Coriolis parameter 
for this case is set to zero, and the formation of hydraulic jumps is to be 
investigated for different values of the mean flow and peak of the 


topographic ridge. 


Il, THE FINITE. ELEMENT METHOD 


A. HISTORICAL INTRODUCTION 

The finite element concept is to a large extent physical rather than 
abstract in nature, and it has been used in a variety of forms for centuries. 
The basic idea has always been to replace an actual problem by a simpler 
one. In other words the finite element method is the replacement of 
continuous functions by piecewise approximations, usually polynomials. 
Indeed the early geometers used "finite elements" to determine an 
approximate value of x. They did this by bounding a quadrant of a circle 
with inscribed and circumscribed polygons, the straight-line segments 
being the finite element approximations to an arc of the circle. In this way 
they were able to obtain extremely accurate estimates. Upper and lower 
bounds were obtained, and by taking an increasing number of elements, 
monotonic convergence to the exact solution would be expected. 
Archimedes used these ideas to determine areas of plane figures and 
volumes of solids, although of course he did not have a precise concept of 
a limiting procedure. The interesting point here is that while many 
problems of applied mathematics are posed in terms of differential 
equations, the finite element solution of such equations utilizes ideas which 
are in fact much older than those used to set up the equations initially. 

The modern use of finite elements really started in the field of 


Structural engineering. Probably the first attempts were by Hrennikoff 


(1941) and McHenry (1943) who developed analogies between actual 
discrete elements and the corresponding portions of a continuous solid. The 
term "finite element" was introduced later by Clough (1960) in a paper 
describing an application in plane elasticity. 

The engineers had put the finite element method on the map as a 
practical technique for solving their elasticity problems, and although a 
rigorous mathematical basis had not been developed, the next few years 
Saw an expansion of the method to solve a large variety of structural 
problems. The workers in the early-1960s soon turned their attention 
towards the solution of non-linear problems. Turner et al. (1960) showed 
how to use an incremental technique to solve geometrically non-linear 
problems, 1. e., problems in which the strains remains small but the 
displacements are large. Stability analysis also comes into this category and 
was discussed by Martin (1965). Plasticity problems, involving non-linear 
material behavior were modelled at this time (Gallagher et al. 1962), and 
the method was also applied to the solution of problems in visco-elasticity 
(Zienkiewicz et al. 1968). 

Finally, besides the static analysis, dynamic problems were also being 
tackled, and Archer (1963) introduced the concept of the consistent mass 
matrix. Both vibration problems (Zienkiewicz et al. 1966) and transient 
problems (Koenig and Davids, 1969) were considered. Thus the period 
from its conception in early-1950s to the mid-1960s saw the method being 


applied extensively by the engineering community. Once it was realized 


that the method could be interpreted in terms of variational methods, the 
mathematicians and engineers were brought together, and many extensions 
of the method to new areas soon followed. In particular it was realized that 
the concept of piecewise polynomial approximation offered a simple and 
efficient procedure for the application of the classical Rayleigh-Ritz 
method. The method had now become an important technique from both a 
practical and theoretical point of view, and the number of published paper 
using this method began to increase at a tremendous rate. 

In the area of meteorology also, the finite-element method has been 
successfully employed in the horizontal representation of atmospheric 
variables in numerical weather prediction and atmospheric modeling 
(Cullen, 1974a, 1974b, 1979; Hinsman, 1975, 1983; Staniforth and 
Mitchell, 1977, 1978). The finite element method when applied to 
meteorological equations gives very accurate phase propagation and also 


handles nonlinearities very well. 


B. FINITE ELEMENT APPROXIMATIONS 

In contrast to the finite difference schemes wherein the domain of 
interest is replaced by a set of discrete points, in the finite element method 
the domain is divided into subdomains called finite elements. The unknown 
function, let us name it u, is represented within each element by an 
interpolating polynomial which is continuous along with its derivatives to a 


specified order within the element. Generally, the interpolating function is 


of lower-order continuity between elements than within an element. Thus 
the fundamental building block in the finite element method is the 


subdomain or element. 


C. THE METHOD OF WEIGHTED RESIDUALS 

There are several ways that can lead us to the same finite element 
formulation. A conceptually simple approach can be formulated using the 
method of weighted residuals. The two primarily special cases of the 
method of weighted residuals (MWR) are the Galerkin and the collocation 
methods. 

In the method of weighted residuals, the desired function u (¢) is 


replaced by a finite series approximation as 


N 
wO=BO)=D/U9,0 (2-1) 
j= 
In general, the set of functions 9; (*), j = 1,2,...,N, can be delineqmages 
both the time and space domain and Uj, j = 1,2...,N, are undetermined 


coefficients. The equation (2-1) can be written 


Neat 
(+) =Uy oy) + Y, UO ©), (2am 


a. 


where 6; (*) satisfy the homogeneous boundary conditions. The functions 
o; (*) are chosen to be polynomials that satisfy certain of the boundary 
conditions imposed on the problem. These functions are variously denoted 
shape functions, basis functions, and interpolation functions, depending 
upon the discipline in which the method is being applied. Even if we 
choose the basis functions to satisfy all boundary conditions, they will not 
normally satisfy the PDE as well. If we now substitute the a (*) into the 


PDE, say Lu - f = 0, we can easily get 
Lu(*)-f=R(e), (2-3) 


where R (*) is a residual. 
Our objective at this point should be to select the undetermined 
coefficients Uj; such that this residual is minimized in some sense. A 


Straightforward scheme would be to set the integral of R (*) to zero as 


[[reavereo (2-4) 


t v 


This scheme, however, generates only one equation for the N unknown 
coefficients Uj. It can be suitably modified by introducing weighting 
mmetions w; (*), 1= 1,2,...,N. 

Setting the integral of each weighted residual to zero yields n 


independent equations 


11 


J [Roe w. (*) dv dt = 0, i= 12 (2455 


t v 


At this point we can solve equation (2-5), in theory at least, for N 
unknown coefficients. Equation (2-5) represents the general equation 
describing the MWR, and a multiplicity of schemes arise out of this one 
expression through the definition of the weighting functions wj. 

Among the MWR family of methods, the Galerkin, subdomain, and 
collocation schemes are most commonly encountered in practice. The one- 
dimensional weighting function for the Galerkin scheme is illustrated in 
Fig. 2.1, for the subdomain scheme in Fig. 2.2, and for the collocation 
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D. GALERKIN METHOD 
The Galerkin method results when the weighting function is chosen to 


be the basis function, as defined in (2-1). Thus we have 


\{ R (*) o, (*) dv dt =0, VS 2 Ne (2-6) 


The basis functions are formally required to be members of a complete set 


of functions. Because a complete set of functions can exactly represent 


2 


basis and weighting 


funcuon 





Fig. 2.1 Schematic representation of the one-dimensional 
weighting functions for the Galerkin method. (It 1s 
assumed that the chapeau function is used as a basis 


fimmetion,):. 
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weighting function 


basis function 


Fig. 2.2 As in Fig. 2.1, except tor the subdomain memoae 
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weighting function 


basis function 





Fig. 2.3 As in Fig. 2.1, except for the collocation method. 
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any function of a given class, the series of (2-1) is inherently capable of 
representing the exact solution as the number of terms in the series 1s 
increased. 

The requirement of completeness allows an alternative interpretation of 
the Galerkin formulation. A continuous function must be zero if it is 
orthogonal to every member of a complete set. The Galerkin method can be 
viewed as a scheme in which the residual is forced to zero in the sense that 
it 1s made orthogonal to the complete set of functions 9; (¢). 

In Fig. 2.1 the weighting function (and therefore the basis function) is 
a hatshaped, piecewise linear function, and because of its hatlike 
appearance, it is sometimes called a “chapeau” function. It is often 
encountered in the formulation of the finite element method. The chapeau 


function is the simplest of the basis functions in common use. 


E. TWO-DIMENSIONAL BASIS FUNCTIONS 

The extension of the weighted residual method to higher dimensions is 
relatively straightforward, provided that regular rectangular subspaces are 
employed. We can easily visualize the two-dimensional case using an 
elementary example, (Lapidus and Pinder, 1982). The discretized domain 1s 
simply illustrated in Fig. 2.4, and the two-dimensional basis function that 


is linear alons each side inte 2-5. 
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Fig. 2.4 Discretized domain for two-dimensional heat flow. 
Finite element nodes are indicated by small circles, and 


the elements themselves by Roman numerals. 
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Two-dimensional basis function that is linear along 


each side. 


Consider now the problem of two-dimensional, time-independent heat 
flow in a rectangular plate with a heat source located at the center of the 


plate. The governing equations for this problem will be given by 


£(T)=T,, + T,,=Q aT) 
T (x,2) = 1 (2S 
T (O,y) = 1 (2-9) 
T, (x,0) = 0 (2-10) 
T, (2,y) =0 (ne) 
Q (xy) = Q, &(x -1) &y - 1), Cay 


where x and y are Cartesian coordinates, Qy is the heat source, and 6 is the 
Dirac delta function. 


ow tep_] 


Let us first define the trial functions: 


9 
T (xy) =t (xy) = Dd) T, 0, (xy). (2-13) 


j=l 
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For convenience, it is advantageous to define basis functions in 


local (€,n) coordinates where upon (2-13) becomes, for gach element, 


4 
T (xy) =f (xy) = > 7,6, Em) (2-14) 


ee 


and 9; (§,1) are bilinear chapeau functions such as those illustrated in Fig. 
aS. 
2,..-Sfep Il 


At this point we can easily formulate the integral equations using 
Galerkin's procedure. We can visualize the whole procedure as the 
requirement of orthogonality between the residual R and each basis 


function, as 


| R (x,y) , (x.y) dx dy = 0, i = 1,...,9, (2-158 
A 


where 
R (xy) =2£°T)-Q. (2 = ice 
A is the domain of integration and £ is defined by (2-7). Substituting (2-7), 


and (2-16) into (2-15) yields 
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J E+ t,- 29, (x,y) dx dy =0, P— eo (2-17) 
A 


Applying Green's theorem to (2-17) to incorporate second- and third-type 
boundary conditions directly into the set of integral equations, equation (2- 


17) becomes 
[164 T, +20) adxdy-[ (TFT, 1, )0,ds=0. 
A S 
fT) (2-18) 


where lx and ly are direction cosines with respect to the normal to the curve 
S, the boundary of the domain A. In this problem, the second term of (2- 
18) will be used to conveniently define the zero-gradient Neumann 
boundary conditions of (2-10) and (2-11). 

If we now substitute the trial functions, as defined by (2-13) into 


(2-18) we can obtain the following set of algebraic equations 


9 
| CX 7, 44: + 5%) + 2, ax dy 


ae y) 'y 
ere 


| (41+ 1,)0, 4s =0. it Nee) (2) 
S 


J 


Because 9; 1s defined such that it is nonzero only over elements adjacent to 
mode 1 (see Fig. 2.4 and 2.5) the integrations of (2-19) may be perionmmen 


piecewise over each element and subsequently summed. Thus we can write 


4 9 
2, J(u (0, 04; + 9; 0,;) + 20, ) dx dy 


e= 


[+t ods=0,  i=1..9, (2-20) 
S 


¢ 


where Ag is the area of element e, and S, is the curve bounding Ae. 


3...-Lfep ITI 
It is now time to formulate the matrix equation. In general, the 
best way is to express (2-20) in terms of the local (€,n) coordinate system 
to facilitate integration (see Fig. 2.6a, Fig. 2.6b). This WS eagime 
accomplished provided that the relationship between derivatives of 9; in 
each coordinate system is readily available. To find this relationship, we 


have to employ the chain rule to obtain 


ae 





Fig. 2.6a Rectangular element in (x,y) coordinates. 
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Fig, 2.6b 
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Same element as Fig. 2.6a transformed into the 
(€,7) local coordinate system. 
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This set of equations can be written as 


(9), a 7 (9), 
oo) Y) (O.), 


OT 


oO) | _,[@ 
EIA ie | (2-21) 
), ,), 


where [J] is the Jacobian matrix. For our problem we can easily evaluate 


the [J] as 
eee 0 

J] = 
0.0 0.5 


and the [J]-! as 
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, [20 0.0 
iar 0.0 2.0] (2 


Equations (2-20) and (2-22) may be combined, which can give us 


4 +1 +] 9 
J J C27, 40,,95, +405 859 #28) =a dn 

-[ (tet) o,ds=0. i = leo ( 22 2a5 
S. 


In changing the limits of integration we introduce the following 


relationship 


dx dy = det [J] d& dn. (2-245 


Because in our problem the elements are all of the same Sizeuuiie 
four integrals appearing in each element in (2-23) are identical fommeas 
element. The assembly process, that is, the transformation from element to 
global integrations, can be easily visualized using Table I which provides 
us with the relationship between Local Node Numbers and Global Node 
Numbers. 

The assembly procedure now calls for extracting information 


concerning the integration at the element level from the Table I) irom 
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Global Nodes 
Element Nodes (Global Node Numbers) 


(Local Node Numbers) | Ef EE 





Table I Relationship between Local Node Numbers and Global 
Node Numbers. 


column I we can see that global nodes 4 and 5 correspond to element nodes 
2 and 3 in element I. In order to find the contribution to the global integral 
concerning these nodes, we have to create Table II which evaluates the 
integrals for one element in local coordinates using equation (2-23). From 
Table II the seeking contribution is found in row 2, column 3, and it is 
equal to -1/6. This value should be placed in matrix location (4,5) of the 
global matrix. However, this value is not final, because there is also 
information to be retrieved from element III. This information is located in 
the position (1,4) of the element matrix of Table II, and it is equals to -1/6. 
This value has to be summed with the previous integral value and the new 
value should to be placed in the same global matrix position (4,5). 
Combination of the two integral values yields the final value of -1/3, which 
can be seen in Table III in row 4, column 5. 

In general, the element coefficient matrix (Table III]) is dittepens 
for each element because of changes in either element geometry aien 
parameter values. Moreover, it is often necessary to perform the 


integrations of (2-23) numerically. 
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Table II Evaluated Integrals for One Element in Local 
Coordinates for equation (2-20). 





Table III Evaluated Integrals for the Four Elements of Figaaee 


in Global Coordinates for Equation (2-20). 
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4... Step IV 


We now solve the matrix equation (2-23) using Table III as 


DemeyG 60.0) 1614 wOs0m 0.0). 0.0 0.0 Tey 
mVicueed/se aie -1/3 -1/3 -1/3 040 0,0 0.0 ee 
OMGMt)/ Gy oee0, 0) -1/3°-1/6 0.0 0.0 0.0 ie 
acme cemonon 92/3 -1/3 1080=1/6 -1/3 0.0 Ty, 
“1/3, -1/3 -1/3 -1/3 8/3 -1/3 -1/3 -1/3 -1/3 | * | Ts 
| 
CuO vee i/o 0.0 =1/3° 4/3 (0.0 -1/3 -1/6 li. 
| 
iMOmOsOMmOs@) 1/6 -1/3 0.0 2/3 -1/6 0.0 ie 
ORCeOMOMNONON 1/3 -1/3 -1/3 21/6 4/3 =1/6 ie 
0.0 0.0 0.0 0.0 -1/3 -1/6 0.0 -1/6 2/3 is 
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Jct etyy 9, ds 
S 
Jct yet o ds 
S 


J ct +t yt yo, ¢s 


S 
0.0 
-Qw 


Jct ety 1,06 4s 
S 


0.0 


™ me ee He a, 


a2 


(2 - Zam 


The line integrals appearing on the right-hand side of (2-25) 


represent flux boundary conditions. In general, when Dirichlet boundaries 


are specified, it 1s necessary to expand and evaluate the line integrals. 


However, when we use functions, that have continuous derivatives up to 


the Oth order as bases, one can condense rows containing the known 


temperature values out of the matrix equation. The reduced matrix equation 


is 


[A] * {b} = 


where 


4/3 
-1/3 
-1/6 


[A] 


-1/3 


(f}, 


-1/3 
8/3 
-1/3 
=1W8: 


(6-26) 
-1/6 -1/3 
WZ 2 
2/3 -1/6|’ aD 
-1/6 4/3 | 
KOE Ob 


33 


Aer O27 So 
Ts5 O25 25 

= ° (2-308 
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Tee a aes 


The coefficients T4, Ts, T7, and Tg represent the values @iueme 
temperature at nodes 4, 5, 7, and 8 because the basis functions are defined 
such that 0; 1s unity at node 1 and zero elsewhere. 


At the same time the analytical solution for our problem is 


T=U+V+W, (2 - sie 
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where 


oo 


2 sin[(n + 1/2)n(@ - y)/a] cosh[(n + 1/2)x(a - x)/a] 


=— SS (7-3) ) 
| a ey (n + 1/2) cosh(n + 1/2) 
: ? £2. y cosh[(n + 1/2)ry / a] sin[(n +1/2)rx / a] (2-33) 
eh (n + 1/2) cosh(n + 1/2) 


Ton 0 (n + 1/2) cosh[(n + 1/2)z] 
(2-34) 


where q is the length of the side of the square (in our case equals to 2), 
and €, n are the x and y locations of the heat source, respectively. The 


analytical solution yields 


qT, 0.756 
T, Cee 
= (2-35) 
JU 0.719 
T, 


0.756, 


a0 


It is obvious that the numerical solution at the singular point (1,1), which 
corresponds to the nodal location 5, is not very accurate. This is not odd 
since we are attempting to represent a rapidly varying function with only 
four bilinear elements. As we mentioned before the series of (2-1) is 
inherently capable of representing the exact solution as the number of terms 


in the series is increased. 
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ROSSBY WAVE 


A. GENERAL 

In order to study motions of atmospheric and oceanic relevance, we 
can use a Shallow, rotating layer of homogeneous fluid which is 
incompressible, and inviscid. This model (Shallow-Water Model) ignores 
completely the presence of stratification, but experience has shown that it 
is capable of describing important aspects of atmospheric and oceanic 
motions. Because of this, it is useful to deal with shallow fluid systems in 
trying to determine general principles of hydrodynamical behavior of the 
atmosphere. 

The major physical characteristics concerning the shallow-water 
model, are found to apply well for more complex systems. If we were to 
analyze the types of wave motions associated with more complex forms of 
the primitive atmospheric equations, (e.g., with vertical stratification, 
compressibility, etc.), we would find essentially the same types, namely 
gravity-inertia waves and Rossby waves, and deep and shallow motions 


would exist simultaneously. 


B. THE SHALLOW-WATER MODEL 
Let us consider a sheet of fluid with constant and uniform density as 


illustrated in Fig. 3.1. The height of the surface of the fluid is given by 


2 


P= CONSTANT 
= 0 





Fig. 3.1 The Shallow-Water Model. 
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h(x,y,t). We also model the body force arising from the potential ® as a 
vector, g, normal to the z = OQ surface. The z-axis coincides with the 
rotation axis of the fluid, so that in this particular case the Coriolis 
parameter f is simply given by 22. The rigid bottom topography is defined 
by the surface z = hg(x,y), which is not a function of time (t). Finally, we 
assume that the fluid is inviscid (u=0), that is, only motions for which 
viscosity is unimportant are considered. 

Also we suppose that a characteristic value for the depth can be 
sensibly chosen, say D, and D also characterizes the vertical scale of the 
motion as well. In the same way we consider that a characteristic horizontal 
length scale for the motion exists, which we call L. The fundamental 


relationship which characterizes the shallow-water theory is given by 


5=- <ay |: (31) 


Here, it is important to note that the fluid is rotating, so that Coriolis 
accelerations can be important. The fluid layer is flat rather than forming a 
spherical shell, and its major physical deficiency as mentioned above is the 
absence of the density stratification which is present in the real 


atmosphere. 
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C. SMALL-AMPLITUDE MOTIONS 
Dealing with the Shallow-Water Model, and since by hypothesis d<<1, 


we are able to express the hydrostatic approximation in the following form: 


Op | 
ee ee (3a) 


We can integrate equation (3.2) resulting in 

p=-pgzt+Ai, y, 0. (3-3) 
Applying now the obvious boundary condition 

p(x,y,h) = Pp, (3-4) 
where po 1S a constant, we can easily get 

p=pg(h-Z)+Ppp. (3. 


At this stage, we can clearly see that the horizontal pressure gradient 


Is Independent of 7, ixc. 


Op | dh 
any ae (oe 
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Senge (3-7) 


in such way that the horizontal accelerations must be independent of z. That 
is why we can also assume that the horizontal velocities themselves remain 
independent of z, if they are so initially. Applying now the Taylor- 
Proudman theorem we are able to write the horizontal momentum equation 


as 


du ou du dh 
3 Ee TV Oy iy = ae (3-8) 
OV OV OV dh 
a Ue (3-9) 


The specification of incompressibility for the Shallow-Water Model, 
decouples the dynamics from the thermodynamics and reduces the equation 
of mass conservation to the condition of incompressibility given in the 


following form: 


a ae 3 - 10 
on a Bylo ra 


Having in mind that u and v are independent of z, equation (3-10) can be 


integrated in z as 
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w(x, y,t) =- (o + ~ + @(x,y,t). (3-Te 


Equation (3-11) can be further manipulated using the condition of no 


normal flow at the rigid surface z = hg resulting in 


dh oO 0 
5 +5, (th hp)u] fg, ee A 2. (3-Ee 


If we now consider the total depth (H) is given by 


7 (3-1 


then the equation of mass conservation (3-12) becomes 


aH 
ot 


+ (UH) += (WH) =0. (3-15) 
Our next step now should be to linearize ate set of equations (3-8), (3- 
9), and (3-15), by studying small-amplitude motions. This is very 
important because the presence of solutions representing free oscillations 
often demonstrates fundamental mechanisms which occur in more 
complicated situations as we mentioned before. 
Let us now consider the thickness of the fluid layer in abSen@euam 


motion to be given as 
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H(x,y,t) = Hy(x% y) +n y.0, (3-16) 


where 7) << Hg. Furthermore 


—— >>u,,. Vu, L217) 


or in other words u and v are considered small enough. Linearizing now the 


equations (3-8), (3-9), and (3-15), we can obtain 


du eo 

me co 7: 
OV _ am 

a = By laid 
on a d _ 

Pero lop o Geo) 


where all the quadratic terms in the dynamical variables u, v, N) with 
respect to the linear terms are ignored. If we now define the linearized 


mass flux vector given by 


U =iU + jv, (3-21) 


where 
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U=uH,, (3-22) 


V=vH,, (3.-238) 


dU 7 on 

Se Se Ole (3-24) 
OV on) 

die wc0N 

ae oe ee (3-2 


At this point it is possible for us to obtain an equation in the single 


variable yj as 


a.a 
= [> +f)n -V. (CV) - gf J (Hymn) =0, (3-2 
aa 
where 
2 
eo (3-28) 
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H dH 
0 on 0 om (329) 


It can be seen that the Jacobian term is the effect of the geostrophic 
wind blowing across isobaths. For low frequency motion and small 
variation in Ho, the first term in (3-27) represents the time rate of change 
of the quasigeostrophic vorticity. The motion that results from this kind of 
balance is a topographic Rossby wave which we examine in more detail in 
the next sections. The velocities components u and v can be also found in 


terms of N, given as 


a an . .om 
fe 8G 5) (3-30) 
a an an 
ta C351) 


D. THE TOPOGRAPHIC ROSSBY WAVE 
Following Pedlosky (1987), let us consider that Ho in equation (3-27) 
varies slightly in the y-direction given by 
H.=D-Ds~ Cree) 
0 9 


L 


where D is a constant, 
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s is the slope of Ho in the y-direction, and 
L is the width of the channel. 

In this particular case (Fig. 3.2), pure geostrophic motion is possible 
only if v is zero. Also, lines of constant Hg are parallel to the x-axis. 
Motion across the isobaths of fluid columns will cause them to stretch or 
contract. Therefore, there is a possibility of a different mode of motion to 
exist depending on the combined effect of rotation and bottom slope. 


Assuming 


g << |, (3-33) 


and 7n to be of the following form 


| = Re { n*(y) expli(kx - ot)] } , (3-345 


substitution in (3-27) yields 


o-F 


gD 


2 Pd 
ain? Sadly * 

-2——+N [ 
ine I dy 





fs 


iy. 
Ss) 
L eG 





-K(1-s2)- k] = 0, (3-250 


with the following boundary conditions: 
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Fig. 3.2 The infinite channel of width L, rotating with angular 


velocity f/2. 
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——+—7 =0, (3-36) 


on y = QO, L. 


If we now assume that 


(l-sZ)=1, (3 - 37m 


which appears to be a very realistic approximation, then (3-35) becomes 


lg * 2 
dy s dn o£ 2 fs : 
_——— -k' -——k = 19) - 
2 eady a C Log da (3 oe 
0 








Here only in terms where s can be compared with quantities of order 


unity has it been neglected. Solving (3-38) we can get 


n = expC) [A sin(ay) + B cos(ay)], (3-30 


where a is given by 


(3-40) 
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Applying now the boundary conditions given by (3-36) we can get the 


following eigenvalue problem 


(0° -f£) (6 - KC.) sin (aL) = 0. (3-41) 


It is obvious, that to the lowest order in s, the slope does not alter the 
Kelvin mode (that is, because the factors multiplying sin (aL) are the same 
as for the case of a flat bottom). The roots corresponding to the zeros of 


sin (aL) are given by 


fksC@ ee 
og AES a2, ney fyi (3-42) 
Lo te 


0 


There are two separate classes of solutions to (3-42) 
a. The first class has frequencies each of which exceeds f. In this 
class the term in s is negligible, and to O(s) we obtain the Poincare modes, 


i S 3 


2 
n v7 


Z 
=~) + O(), neler G3248)) 
L 





o =f +C(k + 


Equation (3-43) shows that the high-frequency Poincare waves are also 


essentially unaffected by the small bottom slope. 
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b. The class having frequencies 6 = O(s), for which the first term 
in (3-42) is negligible, while the second is of O(1). 
The other solution leads us to the dispersion relation for the 
topographic Rossby wave, 1.e., 


fi 
oF =) ’ N= 12555, (3-44) 


ee 
es 


Z 2 
i 1 





obtained by neglecting the first term of (3-41). 


The maximum Rossby-wave frequency should be given by setting 


ma et ole 





k=k =( 5 a) ; (3-45) 
i Cy 
for which we get 
f 
ae max 77 ; (3-46) 
22 £L 1/2 
[In wv + ea 


From equation (3-46) it is clearly seen that the Rossby-wave frequency 
is always less than f. Another important feature of the Rossby wave is that 


its phase speed in the x-direction, which 1s given by 
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sf 

1e 

= 3 - 47 

ae nae ai 
oboe 


Cee 





Z 
0 


is negative for s > Q, and positive for s < Q. Therefore, the wave 
propagates (in the Northern Hemisphere) parallel to the topography, with 
the shallowest water on its right. Also, for high wave number, 1.e., small 
scale, the frequency o decreases with increasing wave number as shown in 


me, 3,3, 


The dynamical fields for the Rossby-wave to lowest order are given by 


N=N, sin) cos(kx - ot + 6) + O(s), (3-48) 
u=- Sn, cos(**) cos(kx - ot + 9) + O68), (3-49) 
v=- Skt sin) sin(kx - ot + 9) + O(5), (3-50) 


where small terms of O(s) have been neglected consistently. 
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Fig. 3.3 A schematic representation of the dispersion diagram 


for Poincare, Kelvin, and topographic Rossby waves 
in a channel: 


From equations (3-48), (3-49), and (3-50), it follows that to lowest 


order in s, and therefore o/f, the fields of motion in the Rossby wave 


Satisfy 
g on 
merere ? 
aa eae ey) 


which are the time-independent forms of (3-31) and (3-32). Because of the 
close relationship between 7, and the pressure field, equations (3-51) and 
(3-52) represent the geostrophic relation for the horizontal motions. 

At this stage, it is clearly seen that to lowest order in 6/f, which 
plays the role of the Rossby number here, the velocity fields, though 
changing with time, remain continuously in geostrophic balance with the 
pressure field. However, the flow is not exactly geostrophic, for then the 
flow would be restricted to travel parallel to the isobaths, i.e., v would 
vanish. Even though the velocity fields are geostrephic to lowest order, it 
is the very small departures from geostrophy that give rise to the wave. It 
is the small cross-isobath flow, which is a nongeostrophic effect, which 
produces the oscillation. The Rossby-wave, whose existence requires both 
s and f to be non zero, is a low-frequency wave oscillation; its period is 


greater than a rotation period. 
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IV, HYDRAULIC JUMPS _IN_ ROTATING 
-AND_NONROTATING SYSTEMS 


A. GENERAL 

A strong and relatively warm wind known as chinook occurs from time 
to time in areas along the eastern slope of the Rocky Mountains. Similar 
phenomena are often observed in the Owens valley on the eastern side of 
the Sierra Nevada. Also cold fronts, approaching the area of Alps, moving 
from north or west, many times undergo severe deformation. This 
deformation may result in a number of important weather events including 
blocking and splitting of air flow on the upstream side or even triggering of 
lee cyclogenesis in the downstream flow. 

Tepper (1952) has proposed that squall lines are modified hydraulic 
jumps, but the idea that there is a certain link between downslope winds 
and hydraulic jumps was proposed by Long (1953). Long suggested that 
the mechanism that produces the hydraulic jump may be similar to one that 
produces strong waves and downslope winds in the atmosphere. Houghton 
and Kasahara (1968) have also proposed other mesoscale hydraulic 
analogies. Williams and Hori (1970) observed a delay in the formation of 
hydraulic jumps in case the Rossby number was less than Q.1 in a rotating 
system. 

However, it has been difficult to confirm this hypothesis about jump 


formation, because there are significant differences between the atmosphere 
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and the simple fluid systems used in the hydraulic theory. One of the major 
reasons for this uncertainty appears to be the vertical propagation of the 
wave energy, which may occur in the real atmosphere but not in fluids 


bounded by a rigid or a free surface. 


B. HYDRAULIC JUMPS IN A NONROTATING SYSTEM 
Let us first consider the one-dimensional Shallow-Water Model. Let us 
also consider a mean flow over an isolated rigid orographic ridge as shown 


in Fig. 4.1. The horizontal momentum equation is given by 


du ou 
+u 


7) 


The continuity equation also, is given by 


—+ua—+h==0, (4-2) 


where h denotes the depth of the fluid, and 
hm is the height of the rigid ridge, which is a function of x. 
At time t = 0, the fluid is set in motion from rest so that for infinite x, 
it has a constant zonal flow uo. After sufficient time has elapsed, the 
solution in the neighborhood of the rigid ridge would be given by the 


steady state solutions of equations (4-1) and (4-2). 
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Fig. 4.1 A cross section view of the one-dimensional Shallow- 


Water Model with a mean flow u. 
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Following Houghton and Kasahara (1968), the steady state solutions for 


the variables u and h are given by 
u + 2gh + 2gh,, =C,, Ae 


CAMA) 


where C,; and Cz are constants. If a flow without a hydraulic jump is 
considered, both C; and C2 are determined by the velocity ug, and the 


height hg of the flow in the region far from the ridge, so that 


z 
C, =u, + 2gh,, (4-5) 


1 


(4x6) 


At this stage we define the following dimensionless parameters Fo, F, 





R, and Us as 
u 
Fj=——, (4-7) 
gh, 
2 
2 Fo 


a7 


__M , 
— (4-9) 
(4-10) 


If we eliminate h from equation (4-3) by using equation (4-4), we are able 


to obtain 


(F=1)U, (RF ce ee (4 oe 


Assuming that the mean flow ug is always positive, then F > 1. Dividing 
equation (4-11) by (F2 - 1) we can get 
2 
R-F 1 
Crees 


Fae ae eal 








Une 270) (4-12) 


Equation (4-12) is a very important relation among the variables Us, 
R, and F. We can easily plot the solution U* to equation (4-12) for given 
values of F2. For instance if F2 equals to 1.02 (Fo = 0.2), equation (4-12) 


becomes 


U2 + (SOR - 51) U, +50 =0. (4-13) 
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We can consider equation (4-13) as determining the values of R that are 
possible with a given value of U+, as shown in Fig. 4.2. It is interesting 
here to note the singular behavior near R = Q. In case of F = 1.125 (Fo = 


0.5), equation (4-12) becomes 


U> + (8R - 9) U, +8 =0, (4-14) 


which leads to Fig. 4.3. Also, if F = 3.0 (Fo = 2.0), equation (4-12) can 


be written in the following simple form 


U> + (0.5R - 1.5) U, + 0.5 =0, (AS) 


which gives Fig. 4.4. 
It 1s very important here to note, that if R has values lower than 


Periieals Where Roritical 1S given by 


R Sipeemnssoo (aay, (4-16) 


critical 


then three real roots exist, as we clearly see in Figs. 4.2, 4.3, and 4.4. In 
iiewcase where R is greater than Re,;jticaj, no physically meaningful solution 
exists (no real roots). In other words, in this particular case we expect the 


occurrence of a hydraulic jump. 
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REPRODIICED AT GOVERNMENT EXPENSE 





Fig. 4.2 Parameter R as a function of solution Us, as given by 
equation (4-12) for F2 = 1.02. 
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Fie. 4.4.As in Pigs 4.2 excep mio etau: 


If we plot equation (4-16), as shown in Fig. 4.5, we are able to find three 
distinct domains. In both domains I and III, the parameter R (function of 
x), has values lower than Reriticay. In this case, a real solution of equation 
(4-16) there exists, which is physically meaningful. On the other hand in 
domain II, the parameter R is greater than Reriticaj, so no physical solution 


exists. 


C. HYDRAULIC JUMPS IN A ROTATING SYSTEM 

Although none of the numerical experiments concerning hydraulic 
jumps (results section), involves rotation, it 1s useful to explore the 
particular effects of rotation on the formation of hydraulic jumps. The basic 
equations for a homogeneous, one-layer, inviscid fluid (Williams and Hori, 


1970), are given by 


du du oh 

ox *ox aa 
Ov OV 

a tua, t a0, (4-18) 
dh oh du 

om on. oe 


where h is the depth of the fluid. 
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4.5 Classification of asymptotic flow conditions as a 
function of the maximum height of the topographic 


ridge, Reritica}, and initial flow speed parameter F. 


We perform a scale analysis by expressing the independent variables 


as 
t=T Tt, (4-20) 


x=Lx’. (4-21) 


The dependent variables are also broken up and scaled as 


(0 Ones (4-22) 

vay ae (4-23) 
h 

See) ——— hy (4-24) 
ne & 


where hm represents the mean depth of the fluid. Also, the Rossby number 


and the Froude number are defined as 


Oh 
= — 4-2 
Ry FL ( 5) 
and 
F A a (4-26) 
12 
(gh_) 
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In order to include the results obtained by the previous analysis which 
are valid for a nonrotating system, we examine the case where Fo < 1 and 


Fo S Ro. The appropriate time scale for doing this is 


i (4-27) 


2 
(gh_) 
and the appropriate v scale 


a 


= 2 - 
V "ia (4-28) 


The nondimensional equations for this case are given by 


2 

du’ du’ oh’ Fo ww 

A aoe ax’ ox ga’ ~” (4-29) 
0 

av av 

oY + Feu’ + u'= 0, (4-30) 

dh’ Poly). 20u (esos 

Fyn Sie 1 a) ae (4-31) 


It is obvious that hydraulic jumps can be formed through the action of the 
nonlinear terms in equations (4-29) and (4-31). Even when Fo is small 


(<<1), they can produce a jump in a nonrotating system. If the Coriolis 
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term in equation (4-29) dominates then a hydraulic jump may be prevented. 


This implies that 


2 
Eo 
— >> F. (4-32) 
R 


The form of the curve dividing the jump region from the nonjump 


region should be given by 


(4-33) 


where A is a constant. The numerical solutions show that the range for A 1s 
from 6.0 to 7.5, which appears to be in agreement with Houghton's 
analytical curve corresponding to A = 6.5. For Fo ~ 1, the above scale 
analysis does not apply, since Fg is then greater than Ro. 

When both Fo ~ I, and Ro ~ 1, all terms in the equations (4-29), (4- 
30), and (4-31) are of the same order. In this case, jumps are expected to 
form. On the other hand, if Ro << 1, the proper time scale is 1/f, while the 
proper v scale should be V = U. The nondimensional equations can be then 
rewritten as 


Rg go yy no, (4-34) 
ot 0 x x 
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t ad ts = 
PTO po (4 338) 
du' OL eee oe aca 
or * Role St Ox * Fax |= i“ 


If we now neglect all terms in Ro, the resulting equations describe an 
inertial oscillation in u and v. We do not expect a hydraulic jump to form in 


this case except perhaps after a very long time. 
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V. MODEL DESCRIPTION 
A. GENERAL 

There are two possible choices of increasing resolution, where desired 
or required, concerning a triangular subdivision. We can use _ near- 
equilateral triangles (Cullen, 1974b) or equilateral triangles (Hinsman, 
1975). Both have the advantage of almost perfect wave propagation 
characteristics. For the same problem we can also use _ rectangular 
subdivisions, although it is not obvious which subdivision is most 
suitable. However, a major advantage of the rectangular subdivision 
(shown in Fig. 5.1) is that it allows algorithms to be developed which take 
full advantage of vector processors. The interesting point here is that the 
Galerkin method has to utilize efficient numerical techniques to be 
considered as a viable option for numerical weather prediction. 

There are several solution procedures available for the Galerkin 
method. No particular attempt is made to optimize computational efficiency 
in this research model. A direct solver is employed using a Gaussian 
elimination procedure. The matrices from the Galerkin procedure are 
decomposed into upper and lower block tri-diagonal matrices. A 
preprocessing, representing the forward substitution stage, can be done 
once. Any time a solution is desired, a back substitution has to be 
performed. That is why the required coefficients for the backward step 


must be stored in a very efficient manner. This particular algorithm 
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Fig. 5.1 Rectangular uniform subdivision for a channel in 


Cartesian coordinates. 
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represents a ‘skyline’ solver, referring to the compact method of storing 
only those coefficients required. This method has both the desired level of 


accuracy and a high degree of computational efficiency. 


B. EQUATION FORMULATION 
In order to integrate the equations governing the free-surface height 
and velocity of an inviscid hydrostatic incompressible fluid we can write 
do 


do 
rerua ac 


s+ o( SoS) = 0 (Ban) 


pmEou or 8 ox” Ca) 
OV Ov Ov db 
Peo oy Oy (5-3) 


where is the geopotential height, 
u is the east/west component of the wind, 
v is the north/south component of the wind, and 
fis the coriolis parameter. 
We now assume that the geopotential height 6, in absence of motion is 


®. Then, in general 


ce 


O(x, yt) =O(x,y) +O (x,y, t), (Sa 


where @ is the mean, and 


o' is a perturbation from the mean. The governing equations now 


can be written 


o',+ BD + (ug’), + (vo), = 0, (5-59 

u.+o' +K,-vQ=0, (5-6) 

V.+o, +K) +uQ=0, (5-7) 

where 

Ou Ov 

a Te (5-8) 
l 2 2 

K= 5 lu +v), (>- 
dv ou 

CF (5-10 


Because of the rapidly moving gravity waves, the stability condition 
for a numerical integration normally requires a much smaller time step than 


for the simple advection equation since ®!/2 >> U. Similar results may be 


q2 


expected with the the more complete equations actually used in numerical 
weather prediction. Since the gravity waves are usually relatively 
unimportant in large-scale weather forecasting, the small time step required 
for computational stability increases the computing time considerably with 
little or no compensation by way of increased accuracy, perhaps even a 
loss. On the other hand, implicit differencing schemes, which may have no 
restriction on the size of the time step, have the serious disadvantage of 
requiring the inversion of a large matrix. 

A semi-implicit scheme has the great advantage of permitting a 
relatively large time step without unduly increasing computation time. In 
other words the semi-implicit scheme slows artificially the propagation 
speed of the fastest gravity waves, which allows a much larger time step 
than the normal Courant-Fredrich-Lewy (CFL) stability criterion and also 
offsets some of the extra computational expense required to solve the 
System of equations assembled at each step. For this reason it is now 
necessary for the implementation of the semi-implicit time discretization to 
rewrite our equations in terms of a velocity potential y and a 


Streamfunction w defined by 


U=X,.-Wy, (5-11) 


Very +Yy, (5-12) 


Ue 


with the following diagnostic relations 


Bai) Ra (5 - ton 
RE eal ie (3 - Tas 


In other words, we can use vorticity and divergence or velocity 
potential and streamfunction as variables instead of velocity components. 
This means that second-order derivatives appear in the equations and 
Poisson equations have to be solved. Using linear elements, the scheme 
obtained for the Poisson equation is very similar to the finite difference 
scheme and can be inverted by the same technique. Cullen and Hall (1979) 
showed that the accuracy of the Galerkin finite element method solution 
was more accurate for the vorticity-divergence formulation of the shallow- 
water equations than for an increase in resolution with the primitive 
formulation. This unstaggered vorticity-divergence together with staggered 
variable formulation gives the best treatment of geostrophic adjustment for 
small-scale features (Williams and Schoenstadt, 1980). 

Following the vorticity-divergence approach and dropping the primes 


for the rest of the chapter, the equations become 


o, + OD + (ud), + (vd), = 0, (5-15) 
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C.+V 0+ (uQ, + (VQ, = 0, Geis) 
D.+V o+V K-(vQ, +(uQ, =0. O52 00) 


where V 2 is the Laplacian operator and € is the relative vorticity given by 


¢ Sn teem (5-18) 


Now we express the velocity as the sum of the rotational and 


irrotational components 


Ra cae (5-19) 
where 

Vi 8 TE (5-20) 

MG Gp Z15) 


and then we are able to rewrite the equations using 


pa. Ss 


Ox dy x, CE) 


Wes 


i YW, (5 - 238 


obtaining the following 


0, + OV x =- (ud), - (V9), (5-24) 
(V *y), = - (uQ, - (VQ, (5-25) 
(V*x), + VO = (vQ, - (WQ, - (KD, - (Ky. (5-26) 


Manipulating now the last equation (5-26) and taking care of the 


bottom topography (assumed to be not a function of time), we can easily 


obtain 
0, + BV “x =-[uO- o5)],-[V - OB), (5 2mm 
(V “y), =- (uQ), - (VQ, , (5 - 28) 
V (x, +0) =(vQ-K,), -(UQ+K,), (5-29) 


where 0g 1s the bottom topography defined by the rigid surface z = hg (x,y) 


not a function of time (t). 
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We now define the domain of integration to be a channel encompassed 
by solid north-south walls with east-west cyclic boundary conditions. The 


boundary condition at the walls should be 
V ~—n= Q, ( 5 aa 30 ) 


where n is the outward pointing normal vector. 

Here it is interesting to note that, in rewriting the equations in this 
form, we have increased their order in x and in y from first to second order 
and should expect that it may be necessary to impose further boundary 
conditions. However, since we have sufficient boundary conditions, 
already, any further specification must not be arbitrary but should be a 
consequence of the previous formulations. Along the walls, the v 


component of the velocity has to be equal to zero, resulting in 


6, =~ fu. (5-31) 


The zonal and meridional components of the wind can be written now as 
USO ay ae (Soe ) 


VEY, +X. (5 - 33) 
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Since the v component of the velocity has to be zero along the north 


and south walls then the obvious boundary condition should be 


a ce (5 - 349) 


and we can satisfy the above condition by simply setting 


y=0, (5 - 35m 


when solving the vorticity equation, and 


N= Os (5-36) 


when solving the divergence equation. As a matter of fact this is an 
Ooverspecification but equation (5-34) would be difficult to apply. Our 
initial conditions of course, must be specifically selected to satisfy both 
equations (5-35) and (5-36). 

As we mentioned before we use a semi-implicit time discretization 
scheme for reasons of computational efficiency. Basically, this semi- 
implicit scheme is simply a modified leapfrog scheme giving a net saving in 
the computational time required to make a forecast for a given time. The 
way in which this is accomplished is to evaluate certain terms implicitly as 


a mean over times (t - At) and (t + At) rather than at time €. 
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Following this approach we evaluate all the terms on the left hand side 


of the equations as an average at times (t - At) and (t + At), and all the right 


hand side at time t. The prognostic equations then become 


d+ OLV x (t- At) + Vy (t+ AD] =- [uO - dp)], - LV (O- dp)],, 


V “W) = - (uQ), - (VQ), 


V Ux, +9 (t+ At) +9 (t- AD] = [(vQ -K], - [4Q+K 


(5-37) 


(eS: ) 


Ge 39) 


If we now solve equation (5-39) for x (t + At) and substitute into equation 


(5-37) we can finally get the following set of equations 


V9" -—@ — = [(vQ -K,], - ((WQ +K,], - 2428 
® (At) ® (At) 


2 
BY. X(t At) l 


—— {[u(-,)], +1v@- 
- a {lu @- Op], + - o)1,), 


V@ +x) =[VQ-K,], -(WQ+K,], 
V (y,) =- (uQ, - (VQ, , 


where 6” is given by 


fo 


(5-40) 


(5-41) 


(5-42) 


0 =[0@-At)-o(t+An ]/2. (5-43) 


At this stage our initial system of three equations in three unknowns 
has been reduced to two Poisson equations and one Helmholtz equation to 
be solved at each time step. Our first step concerning the solution 
procedure should involve solving the > equation (5-40) for a new value of 
oO”. The second step should be then to solve equation (5-41) for (o" + Xx;,) 
and after substitution for x,. At last we solve equation (5-42) for y;. Our 


history variables are >, u, and v and they are updated after each time step 


as 
(t+ At) =20 - 6 (t- Ad), (5-44) 
u (t + At) =2 At Gao - (Wy + u(t - At), (5-45) 
v(t + At) = 220 cay Ri J By sevatherae) (5-46) 
In other words, numerical integration of the three forecast equations 
involves 


a. solve first a Helmholtz equation for 9, 
b. solve a Poisson problem for y, and finally 


c. solve a Poisson problem for x. 
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The space discretization consists of expanding the dependent variables 
in terms of basis functions defined on a variable mesh, and then 
orthogonalizing the error to the basis using the Galerkin procedure 
described in Chapter II. An appropriate approximating function for the 
rectangular subdivision should be a bilinear function (f). In this case the 


forecast set of equations in Galerkin form become 


jes 


) 
ror [iF towo, 6-2, ep}, 


O O l O O 
.f i> [(uQ, f+ 5K, 6) tas] —— Fs (ul - 05) f, + 5000 - 49) f 





Ieee 


(t - at) ff, + [Vx (- a0 gf £,, (5-47) 
® (At) 


d d 
[V7 Hi58-- Joo, 6) 4- foo, 1g, (5-48) 
2,0 * Q OK 
JV? 50,646 Gis [Lg loo,t,- Gogh 


1 


) dK 
- J LF lua, 4+ S58 (5-49) 
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The integral sign here means an area integral over the domain, the j 
subscript refers to Einstein summation for the dependent variables and the 1 
subscript is the ith nodal equation. Following the integration-by-parts 


procedure the final form of the forecast equations is given by 





2 (wa,-k, Dif e, 





7 ce 9; om 
alone i) i $5 el 





Fz 0,6 + GST e + J {hue o0)], GD 
+ (v0 09)], 5, Ff j— — an éf+ | fuc- an, 
+ v(t At) = f- | 486, (5-50) 
[15,3 Dit 5, GH, Hi --] wo, Gye 
JO, ye, (5-51) 


oY. df soat * 7 
(Bo, G26 a ee on ae o, ffi] = 


J J 
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ikea [(vQ), f,- K, Da} e- Aes [u, 6+, SI}, (5-52) 


The line integral along the north and south walls has been dropped 
from the vorticity equation (5-51), since the value of dw/dt is zero on the 
boundaries (w = constant). Also, the line integral along the north and south 
walls has been dropped from the divergence equation (5-52), because the 
value of the normal derivative along the north/south boundaries is zero. As 
we mentioned before the initial conditions should also satisfy the condition 


that the normal derivative of % along the north/south boundaries is zero. 


C. STABILITY ANALYSIS 

Here we analyze the primitive form of the forecast equations and a 
semi-implicit time scheme, since the results will be identical for the 
vorticity-divergence form, (Hinsman, 1983). The one-dimensional 


equations with a mean flow, U, are given by 


du ad du 

oY tty. ae 
OV Ov 

mo t* ai ite 


= +O——=-U—. (5-55) 
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Evaluating the time derivatives with a centered time differencing, and 
averaging the other terms on the left-hand side between time levels (t + At), 


and (t - At), equations (5-53), (5-54), and (5-55) become 


u(x,t + At) - u(x,t - At) A 1 o(x + Ax,t + At) - d(x - Ax,t + At 


[ 


Zt 2 2 Ax 
a o(x + Ax,t - = Oe - Ax,t - ay --U ues + ae wee - oy 
+ f v(x,t), (5-368 


v(x,t + At) - v(x,t - At) =U V+ OD SS aie ) (5-57) 


Zt PLUS. 


o(x,t + At) - O(x,t - At) AG pts + Ax,t + At) - u(x - Ax,t + At) 
2 At Z Zax 


" u(x + Ax,t - At) - u(x - Ax,t - At) o(x + Ax,t) - O(x - Ax,t) 


— }=-U( —— ). (5-58) 
Assuming now a function, F, given as 
F(x,t) = F’ exp[i(kx + at)], (5 - 508 


and substituting into (5-56), (5-57), and (5-58), we can get 
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2 
> 2 - [f+ Ky CO = 0. (5-60) 


At At 

where 
s = sin(@At) + (k') u At, (6-61) 
c =cos(@At), (5-62) 
(k') = sin), (5:4703)) 


The roots of equation (5-60) are given by 
s = 0, (5-64) 
5” = (fAt)’ + c°(kAt) ©. (5-65) 


Requiring @ to be real, the roots of (5-65) yield the following stability 


criterion 


At < ———— (2:00) 
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-YI. INITIAL CONDITIONS 


A. TOPOGRAPHIC ROSSBY WAVE 
Let us first consider the horizontal momentum equations (3-8), and (3- 


9). If we cross differentiate (3-8) with respect to y, and (3-9) with respect 


to x, we obtain 


ou _ ou du au ov du au ov oh (6 1 
dyot dy Ox. dyox dy dy ay? OY 5 Oy Ox’ ; 


ov cl ciee av avo oy + OU, oh ' (6-25 
Ox dt Ox Ox ee Ox Oy Oxdy ox 5 Ox oy” 





If h is eliminated from (6-1), and (6-2), then 


a av au), 8 a a). 8 ay ayy 
ot ox dy Ox dx dy dy ox oy’ 


du du Ov du. du ov 
ae ay ox Woy, oie (6-39 


Using now the following definition: 


av a 
G=0,=30 75, (6-4) 
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where € is the vertical component of the vorticity, equation (6-3) yields 


d¢_ 9, = du = 
att tS -(C + Oi ws me 7) Coss) 


If we now use equation (3-15), equation (6-4) can be written in the 


following form 








dg _ C+f dH 

dt HAH dt tO) 
or 

a Seto 

x | -- EO) (6-7) 


where f is assumed to be constant. 

At this stage, in order to find the proper initial conditions for the 
particular case of the topographic Rossby wave, we can work out a 
simplified theory for Rossby waves (Phillips, 1965). 

Let us uSe cartesian coordinates and simplify the H variation, since H 


is a function of y, as 


H=D(l-sy), (6-8) 


87 


where D is a conStant, and 
s is the slope of H in the y-direction. 


Expanding the (€ + f) /H term as 


Cte ee Sie (C+f)(l+sy) €+f+Csy+fsy 


a 
te ee Se —<=<=———_$_—— ew 


H  D@ecye D D (0 


and considering the term (¢ sy) to be very small, equation (6-7) becomes 


d 
a (G+f sy) =0. | (6 - Tue 


Equation (6-10) can be also written as 


dg dy _ 

Rate O (6-11) 
OT. 

Ss tsv=0 (6-12) 


Assuming small amplitude motion, and no mean flow, equation (6-12) 


becomes 
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oe fsv=0. (ar 13) 
ot 


In the case of a small Rossby number, we can also use the following 


relationships: 
2 
C=V-w Goes) 
and 
mee (Baie 
Ox 


where y = P / (pf). Equation (6-13) then becomes 


G2 OW, 
or 
Cy 2 Oy, 
.  G (6-17) 


mamere > =f s. 


Assuming a solution of the form 
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W = Wy) explip(x - ct)], (6 - lem 


substitution into equation (6-18) yields the following problem for wy) 


2 
d : ; 
(—<- n’) ( isc) explincx - ct)] + B y (iy) explin(x - et)] =0 (6-19) 
dy 
OT 
d° 2 
oY + (u? + Bye 0. (6-20) 
dy’ : 


Nondimensionalizing the domain of integration (as shown in Fig. 6.1), 


our boundary problem for w(y) becomes 





2 
dy B 
+ - (uy, +) wv, =0, (6-21) 
dy 
with 
vy, (0) =0, (6-22) 
y (a) = w,(a), (6-23) 
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Fig. 6.1 Schematic representation of the non-dimensional 


domain of integration. 


)| 


dy,(a)  dy,(a) 


ay ce (6-24) 

for the lower part of the channel,and 

<2 .4)-By=0 (6-25) 
with 

w,(1) = 0, (6 - 26) 

W,(a) = y, (a), (6-27) 

dy, (a) dy, (a) 

WE MGS as 2 (6-28) 


dy dy ° 


for the upper part of the channel. Here 1), and [2 are the corresponding 


wavenumbers. Also 


B,=s, f, (6-29) 


B, =s, f, (6 - 30) 


oe 


where s,; is the slope of H in the y-direction for the lower part, and 
S2 is the slope of H in the y-direction for the upper part. 
There are two distinct cases. In the first case, the lower part of the 


channel controls the phase speed c. The solutions to this case are 


wy, =Asin(Q,y), (one 


for the lower part, and 


W, =B sinh{A, (1 - y)], (6-32) 


for the upper part. For the particular case where a = 0.5 


item See (6-33) 
B 

-Ap=hyt+—, (6-34) 

(Ay) = 2pp - 22, (6-35) 


Furthermore, the coefficients A, and B are to be determined. In the second 
distinct case, the upper part controls the phase speed c. The solutions to 


this case are 


28: 


W, =C sin[A,(1 - y)], (6 - 36) 


for the upper part, and 


W, =D sinh(a.y), (6-37m 


for the lower part. For the particular case where a = 0.5 


B, =-B,, (6 - 38) 

apie ct, (6-39) 
rae 2 

-(Ay) = + 2-22. (6-40) 


Once more, the coefficients © and D are%to be determined: 
Let us now consider the second case, where the upper part controls the 


phase speed c. Boundary conditions (6-23), and (6-27) yield 


sin{X,(1 - a)] C- exp(A, a) - exp(- A, a)] D =0. (6-41) 


In the same way, the boundary conditions (6-24) and (6-28) also yield 
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- A, cos[A(1 sy) |S Mi [exp(, a) + exp(- s a)] D. (6-42) 
Equations (6-41) and (6-42) lead us to the following eigenvalue problem 
A, [exp(A, a) + exp(- A, a)] sin{A,(1 - a)] 


+h, [exp(A, a) - exp(- A, a)] cos[A,(1 - a)] = 0. (6-43) 


Using relationship (6-40), equation (6-43) becomes 


2-1/2 


Me [exp(A, a) + exp(- de a)] sin({(A,)” = cHale | Ce a) 

+ [(A,)° - 205] [exp(A, a) - exp(- A; a)] 

cos({(A,) - 2y3]'" (1 - a)} =0, (6-44) 
where 

L, = = ; (62455) 


Here L is the horizontal length scale, and W represents the width of the 
channel. For our case we choose a channel 30° x 30° longitude by latitude, 


which gives 


aS 


i yee) Os) atk 
W = 4896083.93 m, 
SO 2 equals to 5.4413981 
We are able to solve (6-44) numerically (see Table IV) using Newton's 


method. Using Table IV, we can determine coefficients C, and D for any 


particular mode. For instance if 


dh, = 9.3191, 


ee PA 


then in order to determine the relation between C, and D, the following set 


of equations has to be solved 
0.4912 C - 105.5791 D =0, (6-46) 


4.5787 C - 984.0783 D =0, (6-47) 


resulting in 
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h > h, 


CO Tie 30.8523 





Table IV Numerical solutions of equation (6-44), obtained using 


Newton's method. 


oT 


0.4912 


=705:5791 (6-47) 
or 
Ca; 
D =0.004652 . 


Here C is chosen to be 1. Finally, the solution to the case where the mpiem 


part controls c, should be given by 


yW, =sin[5.2562(1-y)], O5sys1.0, (6 - 48 ) 


for the upper part and 


YW, = 0.004652 [ exp (9.3191 y ) - exp (- 9.3191 y )], 
0.0<y<0.5, (6-49) 
for the lower part. The graph of this solution is shown in Fig. 6.2. 


Following the same approach, with the use of Table IV, we are able to 


obtain 
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lb 


YY 


¥=0.0 | 2 


16 


Fig. 6.2 Graph of the solution wo(y) where the upper part of 
the channel controls c, obtained from equations (6-48), 
and (6-49). 
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y, (upper) = sin [11.1881 (1-y)], O5sy<1.0, (6-500) 


YW, (lower) = - 0.000716 { exp ( 13.5791 y )= exp = 1357 ia 


00 Sy 0S, (6- S18 
for the case where 
A, = 13.5791, 
A, = 11.1881, 
shown in Fig. 6.3, and 
W, (upper) = sin [ 17.3683 (1-y)], OSsy<1.0, (6-52) 


YW, (lower) = 0.0000506 [ exp ( 18.9967 y ) - exp (- 18.9967 y ) J, 


0.0<y<0.5, (6 530) 


for the case where 


dL, = 18.9967, 
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Y 





Y= 0 Z 


Fig. 6.3 As in Fig. 6.2, except for equations (6-50), and (6- 
Sy): 


ror 


d, = 17.3683, 


shown in Fig. 6.4. 
The corresponding solutions for the case where the lower part controls 


the phase Speed c, are givenubelow 
yW, (lower) = sin (5.2562 y ), 0's y = ior (6-545) 


W, (upper) = 0.004652 {exp[ 9.3191(1 - y)] - exp[ - 9.3191(1 - y)]}, 


0.5<y<10, (6-55) 


for the case where 


d, = 9.3191, 


d, = 5.2562, 


shown in Fig. 6.5, 


yw, (lower) = sin (11.1881 y ), 00s ¥= 05, (6-56) 
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¥=-1.0 
iY 


¥=0.0 &. 


iremGndeeAcminerif.gG.), excep fer equations (6-52), and (6- 
See 
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Fig. 6.5 Graph of the solution w,(y) where the lower part of 
the channel controls c, obtained from equations (6-54), 
and (6-55). 
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W, (upper) = -0.000716{exp[13.5791(1 - y)] - exp[ -13.5791(1 - y)]), 


SES SSIS, Ga ae) 
for the case where 
d, = 13.5791, 
A, = 11.1881, 
shown in Fig. 6.6, and 


W, (lower) = sin (17.3683 y), 0.0<y<0.5, (6-58) 


y, (upper) = 0.0000506{exp[18.9967(1 - y)] - exp[ -18.9967(1 - y)]}, 


OS 7100) (6-59) 


for the case where 


1, = 18.9967, 


A, = 17.3683, 
shown in Fig. 6.7. 
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Fig. 6.6 As in Fig. 6.5, except for equations (6-56), andi. 
orth). 


106 


¥=0.5 


C 


0.0 


i a 


Fig. 6.7 As in Fig. 6.5, except for equations (6-58), and (6- 
9): 
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In our experiments we select initial conditions in the form described 
by equations (6-54) and (6-55). That choice appears to be an excellent one, 
providing us with the most proper modes to observe the topographic 
Rossby wave. It is desired that the initial conditions allow a relative 
amount of control for the input parameters as well as satisfing the boundary 
conditions. 

The next logical step is to determine the analytic expression for the 
streamfunction y. Following closely equations (6-54) and (6-55), the exact 


expression for yw, is given by 


nN 





A.A 2 
Vv => [sin yy] [sin 901 Up (F =Ymig) +P (6 - 60) 


for the lower part (0.0 $ y < 0.5), and 





A fe 4 me 
y = 0.004652 = [exp(A, y) - exp( - A; y)] [sin( -- x)] 
® 
“Um > Ymia) * 5s . ( 6618 


for the upper part (O'S Sy Ss ie0)) where 


A = amplitude of perturbation, 


W = width of the channel (4896083.93m), 
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E—"leneth of the channe! (560535 10.7 5m), 


n = wave number, 

Um = mean flow speed, 

Ymid = Middle point of the channel, 
jee 2562. and 

ee = G9 


The first term in the expressions (6-60) and (6-61) represents the 
perturbation part. The second term is the north/south slope necessary to 
Support a mean flow of U,. The last term plays the role of the mean depth 


term. The geopotential height, 06, is related geostrophically to the 


streamfunction, yw, by the following relationship 


) = f, VY; (6- 62 ) 
resulting in 


fen 
Os: : 
>= = sin(a,y) sin(a,x) - fy) U_(y¥-yi.4) + ®, OE o> ) 


for the lower part (0.0 < y < 0.5), and 
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f, A 


9 = 0.004652 





[exp(a,y) - exp(- a,y)] sin(a,x) 


-£) UL, (¥ - Youg) + ®; (6 - 64) 


for the upper part (0.5 <= y = 1-0) 2 wiiene 





rv 
or (6 - 65) 
2m n 
a, = 7 (6-66) 
a, =A, - (606mm 


The u, and v components of velocity can be derived using the 


following geostrophic expressions 


EF By? (6-68 ) 
1 do 
v= oe (6-69) 


resulting in 
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Aa 
u (lower) = les - aa cos(a, y) sin(a,x), Ota 0.5: 





a 
v (lower) = g sin(a, y) Cos(a,x), Ui0'Ss ya 05, 





and 
A : 
u (upper) = U__ - 0.043352 > [exp(a,y) + exp(- a,y)] sin(a,x), 
0.5 sys 1.0, 
Aa, 
Vv (upper) = 5 [exp(ay) - exp(- a,y)] cos(a,x), US Sy a0: 
The initial vorticity can be also derived, using the 

relationship 


Z 
c=Voy, 
resulting in 


A 
Sie > [(a,)” + (a,)'] sin(a,y) sin(a,x), 


fomine lower part (0.0 s y s 0.5), and 


ete 


(6-70) 


con 1) 


[cor 2) 


aes 7S) 


following 


(6-74) 


Cos) 


A | 
¢ = 0.004652 = [(a,)° - (a))°] [exp(ayy) - exp(- a,y)] sin(a,x), (6-76) 
z (ay 


for the upper parti(Ges = 5a) e 

We set the initial divergence equal to zero, assuming the initial fields 
to be almost geostrophic. The initial fields of geopotential, u-component, 
v-component, vorticity, and divergence for the case where the lower part of 
the channel controls the phase speed, c, and no topography effect is 


involved (Equations 6-60 through 6-76), are illustrated in Fig. 6.8. 


B. HYDRAULIC JUMPS 

The stationary theory predicts the formation of jumps in pairs, one 
upstream and one downstream of the rigid ridge (Houghton and Kasahara, 
1968). The classification shown in Fig. 4.5, is strictly valid for non- 
rotating flows. No equivalent theory exists for flows over mountains in a 
rotating system. Houghton (1969) and Williams and Hori (1970) consider 
the transient motion of a shallow water layer on an f-plane without 
mountains starting from an initial velocity disturbance of magnitude U over 
a length L. 

For our case we consider a nonrotating Shallow-Water system (f = OQ). 
It is desired once more, that the initial conditions allow a relative amount 
of control for the input parameters as well as satisfing the boundary 
conditions. The forecast model history-carrying variables are 0, u, and v. 


The analytic expression for the sreamfunction, y, is given by 
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iis 


yw=-Uy, 00<y<wW, (6-79 


where U is the velocity of the flow in the region far from the ridge. Also 


the initial geopotential height, do, 1s set equal to the mean depth, gH. 


Finally, the initial u- and v- components of the velocity are given by 


=e (6-78 ) 


VO) (6-72 
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VIT. EXPERIMENTS AND RESULTS 


Our first experiment involves bottom topography which is composed 
of two regions. These two regions can provide us with either an east-west 
oriented ridge or valley. We consider conditions with no mean flow, the 
objective being to examine how well our model simulates the topographic 
Rossby wave, by comparison with the theoretical phase speed values. 

Our second experiment is to investigate the ability of the same model 
to create hydraulic jumps analogous to that predicted from the analytical 
approach (Chapter IV). In this case we consider a mean flow forced to 
pass over a topographic ridge which extends north-south across the 
channel. In this experiment we will consider several distinct cases 
corresponding to different discrete domains obtained by the theory (Fig. 


4.5). 


A. EXPERIMENT I 

We perform experiment I using the GFEM rectangular model described 
in Chapter VI. The basic difference between the rectangular and triangular 
models is in the approximating polynomials. The rectangular polynomials 
are bilinear while the triangular polynomials are linear. Many integrals 
require evaluation during the integration process. We could use numerical 


quadrature, however, a more efficient method is available through the use 


1B be) 


of natural coordinates. We can accomplish quadrature with no error by 
formula with this method. A description of the natural coordinate method is 
given in Appendix A for the rectangular discretization. 

The diagram for the grids for the rectangular subdivision is shown in 
Fig. 5.1. The domain of integration is 5,653.5 km in the x-directionmane 
4,896.1 km in the y-direction. The model has 12 increments in the x- and 
y- directions, which gives the model 156 degrees of freedom. The Ax is 
471.1 km, and the Ay is 408.0 km. The vale of the Coriolis parameter is 
taken to be 0.00010284, corresponding to 45.0° N latitude. 

The initial conditions for the case of the topographic Rossby wave are 
described in Chapter VI. A small wave perturbation is added to the 
geopotential field which includes the mean height and the required 
north/south slope in case of non zero mean flow. It consists of a wave with 
a wavenumber one, confined primarily to the south part of the channel 
domain. In our case, we choose a mean depth of 1,000 meters, and no 
mean flow. The motion is confined in a channel with cyclic boundary 
conditions as shown in Fig. 7.1. We examine small amplitude wave motion 
encountering different slopes of the bottom topography in the y-direction as 
illustrated in Fig. 7.2, and 7.3. We can visualize the whole setting as if we 
placed a long triangular mountain with its peak centered in the middle of 
the width of the channel. It is obvious that the slope of the lower part of 
the channel, s,;, corresponds to the south slope of the mountain, while the 


slope of the upper part, s2, corresponds to the north slope. No matter what 
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(8) = 0.00010284 





Pig lees chematic representation Of the domain of integration. 
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Fig. 7.2 Schematic representation of the bottom topograpinna em 


the case of triangular mountain. 
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Fig. 7.3 As in Fig. 7.2, except for the case of triangular valley. 
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the case is, the following relationship holds: 


3 Rete (7-1) 


where s; 1s assumed to be greater than or equal to zero for most of 


the cases. 


B. RESULTS I 

We integrate the forecast equations over a time interval of 96 hours, 
and we plot the results every 48 hours. Our first concern is to examine the 
case of no topography. In this particular case I, the peak of the triangular 
mountain is zero, and both slopes, s; and sz, are equal to zero as we can 
see from Table V. The initial conditions are almost purely geostrophic as 
clearly illustrated in Fig. 7.4. The integration produces forecast fields 
almost identical to the initial fields, as shown in Figs. 7.5, and 7.6. That is 
expected, since no topographic or beta effect is involved. If we now 
increase slightly the magnitude both of s; and sz, as shown in Table V for 
case II, the 48 hour integration does not show any significant change in the 
forecast fields, as we can see in Fig. 7.7. However, the 96 hour integration 
yields a very small tendency for a westward shifting valid for all the 
forecast fields, and for the lower part of the channel. 

Our next step, case III, is to increase the peak of the mountain to 


163.1 m. In both the 48, and 96 hour integrations, a tendency for the lower 


120 







a 


Table V Peak (in geopotential meters and in meters), and south 









Slope values for cases I through VII. 
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Fig. 7.4 Initial conditions for experiment I (cases I through 
VII). Contour intervals are 600 m*/s? for geopotential 
height, and 0.2 m/s for u and v. Nodal points are 


denoted by an x. 


STREAMLINES 


E 
ae 
co 
(i) 
2e 
=) 
em 
t— 
FEL, 
OJ 
f= 
© 
(ans 
© 
tJ 
el 


oms’s ocoS@ “400$'S  98c0S'b Secs 6 
IN (SMILBW) ST y-L 


“}0° 


$.a 


4.J 


3.0 
X-AXIS (METERS) 


2.3 


1.0 


BCOS°< = LSS ESS =O > GES 
pt (SyaLIw) SIxy J 


x x x x x x x x 


6e0S*s «= OkOS'M «=A SS ORaOS * Gos 
sole (Suaigu) sixu 1 


X-AXIS (METERS) 


X-AXIS (METERS) 





Fig. 7.5 As in Fig. 7.4, except for case I, and a 48 hour 


integration. 
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Fig. 7.6 As in Fig. 7.5, except for a 96 hour integration. 


124 


CaS PONENT IAL HElemi Si Remit | NES 
TAU = 48 TRU = 48 


$04S96 $0 4896 


7. Se 


x 


x 


6.5.36 


ay) 


% 
/4 


ey SE 


Y-AXIS (METERS) «10° 


35035 4 £28¢ 


3.0 . 
X-AXIS (METERS) 


TEU = 48 
504896 


7 Sue 
7.SC¥ 


a.site 
6.5038 


s 
A x KM K x 


PY 


oatiey) 


6.6987 
SNE 


i See 
Y AXIS (NETORS) » 10° 


Y AXIS EMETERS) 10° 
3.20 6.506 


3-KSS 4.556 


S ae 


SL ; 


f 
70 tes | “\ 
Aes : | 


, . 
’ 

4 ite as 
> - ~ 


VA ‘ ‘ ry4 


SnG 4.0 §.¢ 
X-AXIS (METERS) “id” 





ieee eens in rig. 725, except tor case IT. 


Ne) 


$04S96 


TRU = % 
AXIS (METERS) 


X 


STREAMLINES 


X-AXIS (METERS) 
4.0 5c. 
x10° 


rire 


X-AXIS (METERS) 


t—~ 
a 
Cc) 
-—~+ 
LJ 
a 
aa) 
Ge 
t- 
gee 
LJ 
= 
& 
Ge 
© 
(J 
co 


CIS £  OL5S'D LENS S «= OLOS Seas ¢ 6.0S'¢ O675'O «KOSS OS} OS'S 
im (Sealy Stix I glx (Sualsw) Sty¥yw I 





Fig. 7.8 As in Fig. 7.6, except for case II. 
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part forecast fields for a westward shifting is clearly noticed, as well as a 
small tendency for eastward shifting corresponding to the upper part fields 
(Figs. 7.9, and 7.10). In case IV, which is our last case of positive south 
slope s;, we increase the peak of the mountain to 489.3 m. This value 
corresponds to almost the half of the mean depth value considered for our 
shallow water approximations. In both the 48, and 96 hour integration it 1s 
very evident the significant westward sifting of the forecast fields for the 
the lower part of the channel, and the eastward shifting of the same fields, 
@haracterizing the upper part, as we can see in Figs. 7.11, and 7.12. The 
results obtained here, will be compared quantitatively with the 
corresponding analytical values. 

At this stage we wish to examine cases V, VWI, and VII, all having 
negative south slopes. The above mentioned cases correspond to cases II, 
III, and IV, but with exactly opposite sign slopes. We can visualize the 
whole setting here as if we placed reverse triangular mountains (valleys) of 
different heights with their lowest points centered in the middle of the 
width of the channel, and their peak values always to be given by negative 
values. The 48 hour integration results for case V do not show any 
Significant change in the forecast fields, as we can see in Fig. 7.13, while 
the 96 hour integration yield a very small tendency for a eastward shifting 
valid for all the forecast fields, and for the lower part of the channel (Fig. 
fo). in case VI, the valley is 163.1 m deep in the center. In both the 


48, and 96 hour integrations, a tendency in the lower part forecast fields 
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Fig. 7.11 As in Fig. 7. d3exGept comcaccun = 
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Fig. 7.12 As in Fig. 7.6, 
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Fig. 7.13. As in Fig. 725, exGepmtorn caseue 
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except for case V. 


Fig. 7.14 Asin Fig. 7.6, 
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for a eastward shifting is evident, as well as a small tendency for westward 
shifting corresponding to the upper part fields (Figs. 7.15, and 7.16). The 
same tendency becomes more evident in both the 48, and 96 hour 
integration, 1. e., the significant westward sifting of the forecast fields for 
the the lower part of the channel, and the eastward shifting of the 
corresponding upper fields, as it is clearly shown in Figs. 7.17, and 7.18. 
At this point, we strongly believe that there is a certain link between 
the presence of topography and the observed shifting in all the forecast 
fields, because no shifting at all is been observed in the flat case of no 
topography. That specific link has to be the topographic Rossby wave 
whose existence requires the topographic effect in a rotating system (with 
constant f ). The reason is that the Rossby wave in general, can exist only 
in the presence of an ambient potential-vorticity gradient. In case I, of 
course, no potential-vorticity gradient is present, and that is why we do not 
observe any sign of the Rossby wave. In case II, III, and IV, the positive 
y-direction (shown in Fig. 7.19), 1s also the direction of increasing 
ambient potential-vorticity. If we consider a fluid column initially at rest, 
but later to be displaced in the positive y-direction, then in order to 
conserve its total potential-vorticity, the wave potential-vorticity Cg - Fno, 
has to be decreased. This will balance the excess of the ambient potential 
vorticity in its new place. There are two ways for this to be accomplished. 
a. The fluid will have a tendency to be squeezed, since its new 


position appears to be shallower than the previous one. This will induce the 
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Fig. 7.16 As in Fig. 7.6, except for case VI. 
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Fig. 7.18 As in Fig. 7.6, except for case VII. 
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Direction of the ambient 


potential-vorticity, and 


also positive y-direction. 


Valid for cases II through IV 





Fig. 7.19 The required ambient potential-vorticity gradient. A 
clue for the physical explanation of the topographic 


Rossby wave oscillation. 
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production of negative vorticity due to the vortex-tube compression. 

b. The act of the squeezing can not be complete, since the upoes 
surface is not restricted or bounded. In this case the column will have a 
tendency to ride at least partially up, the slope resulting in a greater value 
of No in its new position than its neighbors. 

It is important here to note that both effects, Cg < 0, and No > 0, act to 
reduce the quantity Co - Fno. Both effects also, will result in a clockwise 
circulation in the fluid around the column. The.clockwise circulation in the 
fluid column C, will force the adjacent column to its right, R, into deeper 
fluid, and the adjacent column to the left, L, to be squeezed into shallower 
fluid, as shown in Fig. 7. 20. The column R, will become the center of a 
counter-clockwise circulation, while column L, will become the center of a 
clockwise circulation. Both contributions from columns L, and R, will 
force the return of column C, toward its original position. This will result 
in an overshoot due to the column C inertia, and the oscillation will 
continue. This is a very simplified view of the phenomenon but it clearly 
shows a very important aspect. The strength of the restoring mechanism 
depends on the vigor of the circulation induced on neighboring fluid 
columns by the displaced column. Following the same approach, in cases 
V, VI, and VII, the positive y-direction is the direction of decreasing 
ambient potential-vorticity, and everything described above applies in the 


opposite sense. 
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at three successive times. Initially collinear and 


isobathwee< IS displaced 
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positioned along an 


move them as shown. The vorticity induced on L 
and R produces a velocity at C with a tendency to 
restore it to its original position. 
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From a purely theoretical point of view, if we recall equation (3-47) we can 
easily find out that for positive values of slope s the phase speed in the x- 
direction is always negative or, referring to our case, westwards. For the 
opposite case of negative values of s, the propagation of the Rossby wave 
in the x-direction is always positive or eastwards. Also, from the same 
equation (3-47), it is obvious that for increasing values of s (in 
magnitude), the propagation phase speed also increases. In other words, if 
we keep increasing the peak of our triangular mountain we should expect 
the presence of the topographic Rossby wave to become more and more 
evident. In order to examine the actual phase speed in more detail, we can 
Fourier analyze the v-component field, and obtain the wave component 


phase speed every 24 hours. The observed phase speed is then given by 


Cc _L@,- >) . 


F (7-2 


20 (t, : t,) 


where L is the length of the channel. These phase speeds, averaged over 96 
hours for cases II through VII, are given in Table VI, and in Figs. 7.21, 
ad 422, 

The analytic phase speed given by (3.41) can be rewritten in the 


following form: 
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Table VI Estimated, estimated corrected by the free surface term, 












estimated corrected by both the free surface and H 
terms, and observed phase speed values (in m/s), for 
cases I] through VII. 
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Figs) 7222 
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C (observed) 
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Comparison of the observed phase speed values, 


with estimated, estimated corrected by the free 
surface term, and estimated corrected by both the 
free surface and H terms, phase speed values (in 


m/s) for cases II through IV (scatter diagram). 
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Fig. 7.22 As in Fig. 7.21, except for cases V through VII. 
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CS (7-35 





where hm is the mountain height, and the term [(t1)* + (A;1)2] is non- 
dimensional. The analytic phase speeds resulting from (7-3) are also given 
in Table VI. 

In general, the analytic phase speeds are all larger in magnitude than 
the corresponding model phase speeds, obtained from (7-2). The reason 
why this happens is that the analytic theory developed here is based on a 
rigid upper lid, but at the same time, the GFEM model used for joun 
forecasts has an upper free surface. If we wish to include the upper free 
surface effect in equation (7-3), the term [(f9)* / gH] must be added to the 


denominator, resulting in 
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where H represents the mean depth value. 
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The value of the new added term is about the half of the denominator 
value, so that it will reduce the analytic phase speed value by around 30%, 
which would bring C and Cr, into more agreement for most of the cases. 
However, the highest mountain peak case, case IV, or the lowest depth 
valley case, case VII, can not be satisfied by only considering this change. 
The above mentioned agreement could be improved for case IV, by using a 
smaller H value into (7-4), since the true average depth in this case is less 
than 1 km. Using the same arguments for case VII, the improvement of the 
desired agreement could be accomplished by using a larger H value into (7- 
4), since the true average depth in this last case is greater than 1 km. 

tables Viandekigs. /.21, 7.22, 7.25, and 7224, compare the observed 
phase speeds with the various theoretical estimates. As expected, each new 
correction improves the agreement with the model phase speed. The phase 
speed for the shallower topography are extremely accurate. The reason for 
the lack of agreement for larger H values is due to the uncertainty for what 
the correct value of H is to use in equation (7-4). Therefore, for higher 
values of the bottom topography the error is larger because using the 
correct value of H in equation (7-4) is most important. Clearly, experiment 
I shows that our numerical model is able to handle the topographic Rossby 


Wave case extremely well. 
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Fig. 7.23 Asin Fig. 7.21, except for bar representation. 
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As in Fig. 7.22, except for bar representation. 
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GC. EXPERIMENT fT! 

We perform experiment II using the same GFEM rectangular model, as 
in experiment I, the only difference is in the resolution of the model 
model now has 24 increments in the x- and y- directions which gives the 
model 576 degrees of freedom. The domain of integration is 727.6 km in 
the x-direction, and 630.1 km in the y-direction. The value of Coniois 
parameter is taken to be zero for all of the cases, corresponding to a 
nonrotating system. 

The initial conditions for experiment II are described in Chapter VI. 


The boundary conditions in x are chosen to be periodic once more, that is 


u (0) =u (L), (7-39 
> (0) =o (L). (7 08 


These boundaries at x = 0, and x = L, are placed sufficiently far from the 
ridge so that the desired asymptotic conditions are well established in the 
vicinity of the ridge before wave motions are able to be fed back into this 
region by the periodic boundary conditions. The height profile of the 


orographic ridge is given by 


i.. = 


M 


H,, sin (—) fone wae 
(7-7) 


O elsewhere, 
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where Z is the width and Hy the height of the mountain, as shown in Fig. 


fe 25. 


D. RESULTS II 

We integrate the forecast equations over a maximum time interval of 
50.6 minutes and we plot the results att = 5.6, 16.9, 28.1, 39.4, and 50.6 
minutes. The five distinct cases we run (I through V), are summarized in 
Tables VII and VIII. Cases III and IV are expected to produce a hydraulic 
jump because they lie in domain II (see Fig. 4.5). In each case the 
equations are integrated from an initial state where u and h are both 
uniform. The u-field for case I is shown in Fig. 7.26. It clearly shows the 
rapid development of a speed maximum over the center area of the mountain 
and slightly on the lee side. A secondary speed maximum also forms over 
the ridge area and it moves upstream with time. The 0-field, shown in Fig. 
7.27, indicates the earlier development of low pressure-field over the lee 
side of the ridge. The high pressure-field over the east part of this pattern 
has an obvious tendency to move upstream with time. Since the main 
feature of the flow, the speed maximum over the center area of the 
mountain, is quite stable we regard this as a no-jump case in agreement 
with the theory presented before in chapter IV. The u-field for case II, a 
case with slightly higher mountain, and stronger mean flow, shown in Fig. 
7.28, iS generally similar to case I except that the perturbation now 


represents a larger fraction of the mean flow. The same arguments appear 
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Fig. 7.25 Schematic representation of the position of bottom 


topography along x-axis, valid for each node per 
horizontal row, for the 


hydraulic 
(experiment IT). 
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Table VII Froude number (Fo), maximum height of the ridge (R), 
parameter F, mean depth (H), mean flow (U), and 


domain (D), for cases I through V. 
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Table VIII Parameter F, maximum height of the ridge (R), 
domain (D), and classification of the asymptotic 
flow conditions, for cases I through V. 
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to be valid for the o-field also, shown in Fig. 7.29. We also regard case II 
aS a no-jump case in agreement with the theory. 

The u-field for case III is shown in Fig. 7.30. In this case the wind 
maximum on the lee side of the mountain, continues to grow with the time. 
At the same time the o-field, shown in Fig. 7.31, indicates that the low 
pressure-field centered on the lee side, changes rapidly to nearly zero. We 
regard case III as a jump case because it is not approaching steady state. 
The resolution of the model is too poor to allow a detailed description of 
the small- scale jump zone. The theory also indicates that this is a jump 
case. Fig. 7.32 contains the u-field for case IV. The o-field for themsamre 
case is given by Fig. 7.33. The behavior of case IV is similar to case III, 
so that, this 1s also a jump case in agreement with the theory. 

The fields of u and 6 for case V, are shown in Figs. 7.34 and 7.35 
respectively. These field patterns are similar to those of cases I and II, the 
only difference is a downstream shifting for both the u- and 9o- fields. The 
growth in the amplitude of the curves appears to be stabilized, so that we 
regard case V as a no-jump Case in agreement once more with the theory. 

In all the investigated cases (I through V), jump formation is indicated 
by both u and o amplitudes which continue to increase with time. In order 
to study the behavior of each one of the jump cases in more detail, better 
resolution will be required as well as a larger domain to reduce the 


boundary effects. 
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VITI. CONCLUSIONS 


In this thesis we have investigated how well a particular shallow-water 
finite element prediction model handles surface topography. In both 
experiments the flow is confined within an east-west channel with periodic 
boundary conditions. 

In the first experiment the bottom is composed of two regions of 
constant and opposite north-south slope with no east-west variation, so that 
the bottom is either an east-west ridge or an east-west valley. In order to 
have the proper initial conditions, the analytic Rossby wave solutions are 
derived. These are obtained by solving the linearized quasi-geostrophic 
equations with the free surface assumed to be rigid. Each solution has a 
Sinusoidal variation with y over one bottom slope and exponential decay 
over the other slope. With the simple Rossby formula the direction of 
propagation is determined by the bottom slope in the region which has the 
largest wave amplitude. These solutions are similar to the trench wave 
solutions obtained by Mysak etal (1979), who used piecewise exponential 
bottom profiles. The numerical solutions with the initial conditions given 
by the linear solutions produced smoothly propagating solutions. The phase 
Speeds were very accurate when the Rossby formula was corrected for free- 
surface effects and mean depth. 

In the second experiment a north-south ridge was placed across the 


channel with sine-squared east-west variation. The Coriolis parameter was 
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set to zero and the initial u and ® were constant. The theory of Houghton 
and Kasahara (1968) for the formation of hydraulic jumps was reviewed. 
The equations were integrated for five initial conditions. In each case a 
speed maximum occurred over or just downstream from the ridge and low 
feiehts were found on the lee side of the ridge. For three of the cases the 
solutions reached an approximate steady state. These agreed with the theory 
which predicted no jumps. The two other cases lead to increasing winds 
and decreasing heights with no steady state. These were jump cases 
according to the theory. In these two last cases the model resolution was 
inadequate to simulate the formation of the hydraulic jumps in detail. 
Finally, the finite element model performed well for two very different 
topographic effects. Further testing 1s required on the jump cases with 
higher resolution. Furthermore, the effect of the  semi-implicit 


discretization on the hydraulic jumps should be determined. 
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APPENDIX 
NUMERICAL QUADRATURE 


A very fast and efficient method is required for the evaluation of the 
integrals obtained by the Galerkin approach. For the rectangular 
subdivision, integration formulas are based on an orthogonal axis 
transformation. In this case, the integrals to be evaluated contain either 
products of the basis functions, products of derivatives of basis functions, 
or a mixture of both. 

An orthogonal axis transformation will allow quadrature formulas for 
rectangles. We are able to transform the rectangle shown in Fig. A.1 by 


using the following 


C=——, (A-1) 








n= (A-2) 


where the values of €, and n at each corner are shown in parentheses. 


Using f; as a basis function, we can express fj as 


pe ee) CR 
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Fig. A.l Orthogonal axis transformation for rectangular 


integration formulas. 
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Derivatives of the basis functon are given by 


Ox a at 

and 
se (A-5) 
y Don 


The interaction coefficients can now be determined for a derivative or 


straight inner product. The straight inner product is given by 
a 
c= ff f,f,dxdy=ab | | f,f,.40an 
pel 


] 
ab 
=> fcisg, purgoa fry n)(1+nn)an 
-l 
=— (2+= =6.6 (242 =1;0,). (A-6) 


Also, the mixed derivative is given by 
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